##########################################################################
# HANGARTNER, LAUDERDALE, SPIRIG (2025) INFERRING INDIVIDUAL PREFERENCES #
##########################################################################

############################################
########### SCRIPT TO RUN MODELS ###########
############################################

rm(list=ls())


#############
# Libraries #
#############

library(rstan)
library(dplyr)
rstan_options(auto_write = TRUE)
options(mc.cores = 1)


# Load functions and code to fit models 
source("00_Models.R")


########
# DATA #
########

# decisions data (not in replication package)
aadecisionsloc <- "Data/AADecisions.csv"
decisions <- read.csv(aadecisionsloc)

# lawyer/paralegal variable
aadataloc <- "Data/aadata.Rda"
load(aadataloc)

Data <- left_join(decisions, DataL)
rm(DataL, decisions)

# judges data
jdataloc <- "Data/jdata.Rda"
load(jdataloc)

Data$Party <- gsub("_\\d{1,2}$", "", Data$J1)





##################
### Run Models ###
##################


## Definitions

MCMCBurn <- 1000
MCMCSample <- 4000
MCMCChains <- 1


## define sample

# Only take 3 judge panel decisions for those submitted in 2007
Data <- Data[Data$panel3==1,]

# Only cases where judge was at court when case was submitted (in 2007)
Data <- Data[Data$J1 %in% JData$TAG & Data$J2 %in% JData$TAG & Data$J3 %in% JData$TAG,]

# Drop cases where the outcome is NA
Data <- Data[!is.na(Data$granted),]

# Prep Stan Data
JudgeName <- sort(JData$TAG)
PartyName <- sort(unique(JData$Party))
Y <- Data$granted # 1 if Outcome=="Gutheissung" (0 if outcome %in% c("Abweisung", "NE"))
M <- length(Y)
n_judges <- length(JudgeName)
n_parties <- length(PartyName)
Cj <- as.factor(Data$Land)
Cj <- relevel(Cj,"LK") # sri lanka as reference country
Kc <- nlevels(Cj)

Ij <- matrix(NA,M,3) # set up matrix with chair, second, third judge
for (j in 1:M){
  Ij[j,] <- c(as.character(Data$J1[j]),
              as.character(Data$J2[j]),
              as.character(Data$J3[j]))
  } 
Ij.m <- data.frame(factor(Ij[,1],levels= JudgeName),
                   factor(Ij[,2],levels= JudgeName),
                   factor(Ij[,3],levels= JudgeName))
names(Ij.m) <- c("J1","J2","J3")

Pj <- apply(X=Ij, 2, FUN=function(x) gsub("_\\d{1,2}", "", x)) # for party FEs 
Pj.m <- data.frame(factor(Pj[,1],levels= PartyName),
                   factor(Pj[,2],levels= PartyName),
                   factor(Pj[,3],levels= PartyName))
names(Pj.m) <- c("P1","P2","P3")

      
nrow(Data) # 1739 for main specification


###################
# Fit Main Models #
###################


standata <- list(
  Y=Y,
  Ij=cbind(as.numeric(Ij.m[,1]),as.numeric(Ij.m[,2]),as.numeric(Ij.m[,3])),
  Lj=as.numeric(Cj),
  M=M,
  N=n_judges,
  K=Kc
)

        
# chair model
posterior_chair = stan(model_code = stancode_covariate_chair, data = standata, init="random",
                                       seed=2024, pars=c("theta","pi","phi","ll","inconsist"), 
                                       iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, chains = MCMCChains,
                                       refresh=1000, control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))

# median model
posterior_median = stan(model_code = stancode_covariate_median, data = standata,
                                init="random", seed=2024,
                                pars=c("theta","pi","phi","ll","inconsist"), 
                                iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                                chains = MCMCChains, refresh=1000, 
                                control = list(adapt_delta = 0.999, stepsize = 0.01, max_treedepth = 15))


# mixture model
posterior_mixture = stan(model_code = stancode_covariate_mixture, data = standata,
                                     init="random", seed=2024, 
                                     pars=c( "theta","pi","phi","ll","inconsist","inconsist_chair","inconsist_median", "delta"),
                                     iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                                     chains = MCMCChains, refresh=1000,
                                     control = list(adapt_delta = 0.999, stepsize = 0.01, max_treedepth = 15))


#########################
# Fit Additional Models #
#########################

# second
posterior_second = stan(model_code = stancode_covariate_second, data = standata,
                                init="random", seed=2024, 
                                pars=c("theta","pi","phi","ll"),  
                                iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                                chains = MCMCChains, refresh=1000, 
                                control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))

# third
posterior_third= stan(model_code = stancode_covariate_third, data = standata,
                              init="random",seed=2024,
                              pars=c("theta","pi","phi","ll"), 
                              iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                              chains = MCMCChains,refresh=1000, 
                              control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))


# max
posterior_max = stan(model_code = stancode_covariate_max, data = standata,init="random",
                             seed=2024, pars=c("theta","pi","phi","ll"), 
                             iter = MCMCSample + MCMCBurn,  warmup= MCMCBurn, chains = MCMCChains,
                             refresh=1000, control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))

# min
posterior_min = stan(model_code = stancode_covariate_min, data = standata,init="random",
                             seed=2024, pars=c("theta","pi","phi","ll"), 
                             iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                             chains = MCMCChains,refresh=1000, 
                             control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))

# null
posterior_null = stan(model_code = stancode_covariate_null, data = standata,init="random",
                              seed=2024, pars=c("theta","pi","phi","ll"),
                              iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                              chains = MCMCChains, refresh=1000, 
                              control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))
          
              
# mean
posterior_mean = stan(model_code = stancode_covariate_mean,
                                            data = standata, init = "random",
                                            seed=2024, pars = c("theta", "pi", "phi", "ll"), #, "misclass"
                                            iter = MCMCSample + MCMCBurn, warmup=MCMCBurn, chains = MCMCChains, refresh = 1000, control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))
              

# mixture model with party instead of judge FEs
standata_p <- list(
            Y=Y,
            Ij=cbind(as.numeric(Pj.m[,1]),as.numeric(Pj.m[,2]),as.numeric(Pj.m[,3])),
            Lj=as.numeric(Cj),
            M=M,
            N=n_parties,
            K=Kc
          )
          

posterior_mixture_p = stan(model_code = stancode_covariate_mixture, data = standata_p,
                                           init="random", seed=2024, 
                                           pars=c("theta","delta","pi","phi","ll"), #,"inconsist"
                                           iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                                           chains = MCMCChains, refresh=1000,
                                           control = list(adapt_delta = 0.999, stepsize = 0.01, max_treedepth = 15))

          
# chair model with party instead of judge FEs
posterior_chair_p = stan(model_code = stancode_covariate_chair, data = standata_p,init="random",
                                         seed=2024, pars=c("theta","pi","phi","ll","inconsist"), 
                                         iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, chains = MCMCChains,
                                         refresh=1000, control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))
          

# median model with party instead of judge FEs
posterior_median_p = stan(model_code = stancode_covariate_median, data = standata_p,
                                          init="random", seed=2024,
                                          pars=c("theta","pi","phi","ll","inconsist"), 
                                          iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                                          chains = MCMCChains, refresh=1000, 
                                          control = list(adapt_delta = 0.999, stepsize = 0.01, max_treedepth = 15))


# mixture model with sample split by chair judge seniority

judge_seniority <-  cut(JData$seniority, breaks = quantile(JData$seniority, probs = seq(0,1,0.5)), include.lowest = TRUE, labels = c(0, 1))
Senior <- as.character(judge_seniority[match(Data$J1, JData$TAG)]) 

print(paste0("Average years of service Senior Chair Judges: ", round(mean(JData$seniority[judge_seniority==1]),1)))
print(paste0("Average years of service Junior Chair Judges: ", round(mean(JData$seniority[judge_seniority==0]),1)))

standata_s <- list(
  Y=Y,
  Ij=cbind(as.numeric(Ij.m[,1]),as.numeric(Ij.m[,2]),as.numeric(Ij.m[,3])),
  Lj=as.numeric(Cj),
  M=M,
  N=n_judges,
  K=Kc,
  SC=as.numeric(Senior)
)

posterior_mixture_senior = stan(model_code = stancode_covariate_mixture_split, 
                                        data = standata_s,init="random",
                                        seed=2024, pars=c("theta","delta1", "delta0", "pi","phi","ll", "inconsist"), 
                                        iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, chains = MCMCChains,
                                        refresh=1000, control = list(adapt_delta = 0.99, stepsize = 0.01, max_treedepth = 15))


# independent vote mixture model
posterior_mixture_independent = stan(model_code = stancode_covariate_mixture_independent, data = standata,
                                               init="random", seed=2024,
                                               pars=c("theta","pi","phi","ll","inconsist","delta"),
                                               iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                                               chains = MCMCChains, refresh=1000, 
                                               control = list(adapt_delta = 0.999, stepsize = 0.01, max_treedepth = 15)) 


# independent vote chair model
posterior_chair_independent = stan(model_code = stancode_covariate_chair_independent, data = standata,
                                             init="random", seed=2024,
                                             pars=c("theta","pi","phi","ll","inconsist"), 
                                             iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                                             chains = MCMCChains, refresh=1000, 
                                             control = list(adapt_delta = 0.999, stepsize = 0.01, max_treedepth = 15))          


# independent vote median model
posterior_median_independent = stan(model_code = stancode_covariate_median_independent, data = standata,
                                              init="random", seed=2024,
                                              pars=c("theta","pi","phi","ll","inconsist"), 
                                              iter = MCMCSample + MCMCBurn, warmup= MCMCBurn, 
                                              chains = MCMCChains, refresh=1000, 
                                              control = list(adapt_delta = 0.999, stepsize = 0.01, max_treedepth = 15))



posteriorfilelocation <- paste0("SavedPosterior.Rdata")
          
save(list = ls(.GlobalEnv, pattern = "posterior_|^Data$|JudgeName$|PartyName$"), file=posteriorfilelocation)



